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Abstract.  Ising  models,  and  the  physical  systems  described  by  them,  play  a 
central  role  in  generating  entangled  states  for  use  in  quantum  metrology  and 
quantum  information.  In  particular,  ultracold  atomic  gases,  trapped  ion  systems, 
and  Rydberg  atoms  realize  long-ranged  Ising  models,  which  even  in  the  absence 
of  a  transverse  field  can  give  rise  to  highly  non-classical  dynamics  and  long-range 
quantum  correlations.  In  the  first  part  of  this  paper,  we  present  a  detailed  theoretical 
framework  for  studying  the  dynamics  of  such  systems  driven  (at  time  t  =  0)  into 
arbitrary  unentangled  non-equilibrium  states,  thus  greatly  extending  and  unifying 
the  work  of  Ref.  [1],  Specifically,  we  derive  exact  expressions  for  closed-time- 
path  ordered  correlation  functions,  and  use  these  to  study  experimentally  relevant 
observables,  e.g.  Bloch  vector  and  spin-squeezing  dynamics.  In  the  second  part,  these 
correlation  functions  are  then  used  to  derive  closed-form  expressions  for  the  dynamics 
of  arbitrary  spin-spin  correlation  functions  in  the  presence  of  both  T\  (spontaneous 
spin  relaxation/excitation)  and  Ti  (dephasing)  type  decoherence  processes.  Even 
though  the  decoherence  is  local,  our  solution  reveals  that  the  competition  between 
Ising  dynamics  and  T\  decoherence  gives  rise  to  an  emergent  non-local  dephasing  effect, 
thereby  drastically  amplifying  the  degradation  of  quantum  correlations.  In  addition 
to  identifying  the  mechanism  of  this  deleterious  effect,  our  solution  points  toward  a 
scheme  to  eliminate  it  via  measurement-based  coherent  feedback. 
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1.  Introduction 

Interacting  spin  models  provide  a  remarkably  accurate  description  of  a  diverse  set  of 
physical  systems,  ranging  from  quantum  magnetic  materials  [2,  3,  4]  to  quantum  dots  [5], 
nitrogen  vacancy  centers  [6],  superconducting  qubit  arrays  [7],  ultracold  atomic  gases  [8], 
and  trapped  ions  [9,  10].  Despite  being  relatively  simple,  and  often  admitting  accurate 
theoretical  descriptions,  they  support  a  variety  of  complex  equilibrium  properties  found 
in  real  materials,  e.g.  emergent  spatial  ordering  [2],  quantum  criticality  [11],  and 
nontrivial  topological  phases  [12,  13].  While  the  equilibrium  physics  of  the  simplest 
quantum  spin  models  is,  with  many  notable  exceptions,  fairly  well  understood,  the  study 
of  driven,  dissipative,  and  otherwise  non-equilibrium  behavior  is  comparatively  full  of 
open  questions:  Under  what  circumstances  can  equilibrium  correlations  survive  coupling 
to  a  noisy  environment  [14,  15,  16]?  To  what  extent  do  the  concepts  of  criticality  and 
universality  extend  to  dynamics  and  non-equilibrium  steady-states  [17,  18,  19,  20]? 
Can  interesting  quantum-correlated  states  be  stabilized  by  (rather  than  degraded  by) 
decoherence  [21,  22,  23,  24,  25,  26,  27,  28,  29]?  In  recent  years,  it  has  become  increasingly 
apparent  that  non-equilibrium  dynamics  is  ideally  suited  to  investigation  by  quantum 
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simulation  [30],  making  such  questions  especially  timely  and  important.  Moreover,  there 
are  many  examples  where  interesting  non-equilibrium  states  of  matter  are  more  readily 
achievable  than  low  temperature  equilibrium  states  in  ultracold  neutral  gases  [31],  polar 
molecules  [32],  and  trapped  ions  [33]. 

With  these  motivations  in  mind,  in  this  manuscript  we  develop  a  general  formalism 
for  calculating  unequal-time  correlation  functions  of  arbitrary-range  Ising  models  driven 
far  out  of  equilibrium  at  time  t  —  0,  thus  establishing  a  comprehensive  toolbox  for  the 
description  of  non-equilibrium  dynamics  in  a  simple  context.  In  addition  to  providing 
a  tractable  example  of  quantum  many-body  spin  dynamics,  the  Ising  model  is  realized 
to  a  good  approximation  in  a  variety  of  experimentally  relevant  systems.  And,  despite 
its  simplicity,  Ising  spin  dynamics  is  known  to  be  useful  for  the  production  of  entangled 
states  with  applications  in  quantum  information  and  precision  metrology  [34],  Our 
results  constitute  a  unified  approach  to  describing  experiments  aimed  at  producing 
such  states  [35,  36,  37,  38,  39,  33],  and  facilitate  a  quantitative  treatment  of  a  variety 
of  unavoidable  experimental  complications,  e.g.  long-range  (but  not  infinite-range) 
interactions  and  initial-state  imperfections. 

Ultracold  atomic  systems  are  also  well  suited  to  the  controlled  inclusion  of 
dissipation,  prompting  a  number  of  theoretical  proposals  to  exploit  dissipation  for  the 
creation  of  interesting  quantum  states  [40,  19,  21,  22] — remarkably,  such  ideas  are 
already  coming  to  experimental  fruition  [24,  41],  Verifying  that  these  experimental 
systems  behave  in  the  expected  manner  in  the  presence  of  dissipation,  however,  is 
extremely  challenging,  in  large  part  clue  to  the  numerical  complexity  of  simulating 
dynamics  in  open  quantum  systems  and  the  scarcity  of  exact  solutions.  The  Ising 
model,  especially  as  implemented  in  trapped  ion  experiments  [39,  42,  33],  poses  a  unique 
opportunity  to  study  the  effects  of  dissipation  in  a  controlled  and,  as  we  will  show, 
theoretically  tractable  setting.  In  the  absence  of  dissipation,  an  important  issue  in  the 
Ising  model  is  whether  ground  state  correlations  survive  the  application  of  an  equilibrium 
coherent  drive  that  does  not  commute  with  the  interactions — i.e.  a  transverse  field.  In 
the  dissipative  Ising  model  an  analogous  question  can  be  posed:  How  does  the  system 
respond  to  being  driven  incoherently  by  processes  that  do  not  commute  with  the  Ising 
interaction?  Our  formalism  for  the  calculation  of  unequal-time  correlation  functions 
allows  us  to  definitively  answer  this  question. 

Quite  surprisingly,  non-equilibrium  dynamics  in  the  Ising  model  remains  solvable 
in  the  presence  of  non-commuting  dissipation  [1],  even  for  completely  arbitrary  spatial 
dependence  of  the  Ising  couplings  (and  therefore  in  any  dimension).  This  manuscript 
substantially  extends  the  groundwork  laid  in  Ref.  [1],  where  the  quantum  trajectories 
technique  was  used  to  obtain  a  closed-form  solution  of  non-equilibrium  dynamics  for  a 
special  class  of  initial  states.  The  present  work  not  only  provides  a  more  direct,  unified, 
and  comprehensive  exposition  of  the  relevant  theory,  but  also  generalizes  those  results 
to  include  a  broader  class  of  initial  conditions  and  observables,  and  applies  the  solutions 
to  a  number  of  experimentally  relevant  problems  (most  of  which  had  previously  been 
explored  only  by  numerical  or  approximate  techniques,  if  at  all).  We  also  develop  a 
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clear  physical  picture  of  the  interplay  between  coherent  interactions  and  spontaneous 
spin  flips,  which  reveals  that  Tf  decoherence  is  much  more  detrimental  to  entanglement 
generation  than  might  be  naively  expected.  However,  our  solution  also  points  toward  a 
measurement-based  feedback  scheme  that  can  mitigate  its  detrimental  effects. 

The  organization  of  the  manuscript  is  as  follows.  In  Sec.  2  we  consider  the  coherent 
(Hamiltonian)  far-from-equilibrium  dynamics  of  an  Ising  model  with  arbitrary  spin- 
spin  couplings.  Our  results  comprise  a  unified  framework  for  calculating  unequal-time 
correlation  functions  starting  from  arbitrary  unentangled  pure  states.  As  special  cases, 
these  results  will  be  applied  to  calculating  Bloch-vector  dynamics,  arbitrary  equal  time 
spin-spin  correlation  functions,  and  two-time  dynamical  correlation  functions.  These 
results  are  substantially  more  general  than  any  already  available  in  the  literature 
[43,  44,  34,  45,  1],  and  will  help  quantify  the  quantum-enhanced  precision  in  metrology 
experiments  using  trapped  ions  and  ultracold  neutral  atomic  gases.  In  Sec.  3 
we  consider  the  effect  of  Markovian  decoherence  on  this  dynamics,  incorporating 
dephasing,  spontaneous  excitation,  and  spontaneous  relaxation.  Because  the  excitation 
and  relaxation  processes  do  not  commute  with  the  Ising  dynamics,  including  them  is 
especially  nontrivial:  we  work  in  the  interaction  picture  of  the  Ising  Hamiltonian  and 
incorporate  them  as  time-dependent  perturbations.  Terms  in  the  perturbative  expansion 
are  evaluated  using  the  tools  laid  out  in  Sec.  2,  and  by  summing  the  perturbation  theory 
to  all  orders  we  obtain  an  exact  (closed-form)  description  of  the  dissipative  dynamics 
of  arbitrary  two-point  correlation  functions.  This  is  in  stark  contrast  to  the  behavior 
of  a  coherently  driven  Ising  model,  where  such  a  perturbative  expansion  cannot  in 
general  be  resummed.  An  interesting  feature  revealed  by  our  exact  solution  is  that  the 
spin  dynamics  undergoes  an  oscillatory-to-damped  transition  at  a  critical  dissipation 
strength,  which — in  the  absence  of  a  coherent  drive — cannot  occur  at  the  single-particle 
or  mean-field  level.  This  feature,  along  with  the  more  general  structure  of  our  solution, 
is  demonstrated  by  solving  for  spin  dynamics  in  a  nearest-neighbor  Ising  model.  We 
conclude  this  section  by  casting  our  solution  in  terms  of  a  clear  physical  picture,  in 
which  Ti  decoherence  (spontaneous  excitation/relaxation),  through  its  interplay  with 
the  Ising  dynamics,  gives  rise  to  an  emergent  non-local  dephasing  process.  Section  4 
applies  the  solution  to  calculating  experimentally  relevant  observables  in  a  dissipative 
version  of  the  one-axis  twisting  model.  We  show  that  the  emergent  dephasing  discussed 
in  Sec.  3  severely  diminishes  the  precision  enhancement  achievable  compared  to  that 
obtained  in  the  absence  of  decoherence.  However,  we  also  show  that,  in  special  cases, 
this  degradation  can  be  prevented  by  a  measurement-based  feedback  mechanism.  In  Sec. 
5  we  summarize  our  results  and  pose  a  number  of  unanswered  questions  that  would  be 
interesting  to  address  in  future  work. 

2.  Coherent  dynamics  in  Ising  models 

Our  goal  is  to  develop  a  unified  strategy  for  describing  the  dynamics  of  a  collection  of 
spin-1/2  particles  interacting  via  Ising  couplings  and  initially  (at  time  t  =  0)  driven  far 


CONTENTS 


5 


out  of  equilibrium.  In  the  absence  of  a  magnetic  field,  the  most  general  form  for  an 
Ising  model  isf 

«  =  E  a) 

j<k 

where  crj  are  z  Pauli  matrices  and  the  indices  j,  k  label  lattice  sites  located  at  spatial 
positions  ry.  The  coupling  constants  Jjk  are  left  completely  arbitrary,  and  hence  there 
is  not  necessarily  any  notion  of  dimensionality.  In  many  physical  realizations  of  this 
Hamiltonian,  such  as  trapped  ions,  neutral  atoms,  or  Rydberg  atoms,  the  couplings 
exhibit  a  roughly  power-law  spatial  dependence,  Jjk  =  J\rj  —  Vk\~^- 

Because  there  is  no  transverse  field  (a  term  oc  h  cr)),  the  eigenstates  of  TL  can 
always  be  chosen  to  be  simultaneous  eigenstates  of  all  the  <rj  (with  eigenvalues  a J  =  ±1). 
As  a  result,  the  partition  function  Z(j3)  =  tr  [e_/W]  (with  the  inverse  temperature) 
is  identical  to  that  of  a  classical  Ising  model,  and  the  equivalence  of  all  equilibrium 
properties  follows.  This  is  the  sense  in  which  the  Ising  model  without  a  transverse 
held  is  often  said  to  be  “classical”  (even  though  it  is  a  quantum  Hamiltonian  acting 
on  vectors  in  a  Hilbert  space).  In  passing  we  note  that  classical  Ising  models  can,  of 
course,  be  highly  nontrivial:  for  example,  disordered  or  frustrated  couplings  give  rise  to 
classical  glassiness  [46,  47,  48].  Out  of  equilibrium,  however,  this  notion  of  classicality 
is  inapplicable.  While  a  thermal  density  matrix  p(/3)  =  Z~1e~l3n  commutes  with  TL, 
the  density  matrix  describing  some  non-equilibrium  initial  conditions  will  not  in  general 
commute  with  the  Hamiltonian,  and  nontrivial  dynamics  will  ensue.  This  dynamics — 
which  has  no  direct  analogue  in  the  classical  Ising  model — is  generically  characterized  by 
the  growth  of  entanglement,  leading  in  some  special  cases  to  spin-states  with  applications 
in  quantum  information  and  precision  metrology [34], 

Everywhere  in  this  manuscript,  we  assume  the  system  starts  in  a  pure  state  that 
is  a  direct  product  between  the  various  spins  [Fig.  1(a)]. §  The  most  general  such  state 
can  be  specified  by  choosing  spherical  angles  9j  and  <J)j  describing  the  orientation  of  the 
spin  at  each  site  j  [Fig.  1(b)].  Defining 

fi(  1)  =  e“^/2  cos  |,  fj(- 1)  =  e^/2  sin  (2) 

and  states  \<jj)  that  are  eigenstates  of  <rj  with  eigenvalues  Gj  =  ±1,  such  a  state  can  be 
written 

U(*  =  0))=  (3) 

=  (4) 

3 

For  uniform  d3  =  6  and  4>j  =  (j),  the  state  |T(0))  is  frequently  encountered  in  experiments 

t  Note  that  many  references  studying  long-ranged  Ising  models — i.e.  those  for  which  JT  Jj  is  an 
extensive  quantity — often  normalize  the  interaction  by  dividing  by  the  number  of  spins  N .  We  drop 
this  constant  here  to  avoid  cluttering  the  notation. 

§  Unentangled  but  mixed  initial  density  matrices  can  be  easily  accounted  for  by  suitable  averaging  of 
the  expressions  given  over  the  initial  ensemble. 
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Figure  1.  The  only  restriction  on  the  initial  state  is  that,  it  be  unentangled  (i.e. 
a  product  state  over  the  various  sites  in  the  spin  model).  For  example,  one  can 
imagine  (a)  an  initial  state  with  some  slow  spatial  variations  in  the  spin  angles  due 
to  inhomogeneities  in  the  pulse  strength  of  a  Ramsey-type  experiment.  The  notation 
used  to  characterize  the  state  of  any  one  spin  is  shown  on  the  Bloch  sphere  in  (b). 


implementing  Ramsey  spectroscopy  [49,  33,  32],  and  spatially  varying  angles  could  be 
used,  for  example,  to  describe  the  effects  of  defects  or  excitation  inhomogeneities  in  such 
experiments  [50,  51,  52], 

Essentially  all  properties  of  the  non-equilibrium  dynamics  axe  contained  in  unequal- 
time  correlation  functions  of  the  spin  operators  af  and  a~  (these  subsume,  of  course, 
the  time  evolution  of  all  equal-time  correlation  functions).  We  focus  first  on  the  case 
where  only  operators  af  occur 


S  =  { Tc  («£(«)•■•  (O  ■  •  ■  <(fi))>  • 


(S) 


Here  a,  b  =  ±,  and  the  time  dependence  of  the  operators  is  given  in  the  Heisenberg 
picture  of  H 


=  eimcra-e~itn .  (6) 

The  time-ordering  operator  7c  orders  all  operators  along  a  closed-time-path  C  shown 
in  Fig.  2,  with  times  t  occurring  on  the  forward  path  and  times  t*  occurring  along  the 
backwards  path.  This  closed-time  path  ordering  occurs  naturally,  for  instance,  in  any 
perturbative  treatment  of  additional  non-commuting  terms  in  the  Hamiltonian.  In  Sec. 
3  we  encounter  this  situation  when  treating  a  dissipative  coupling  to  an  environment  , 
but  the  same  structure  occurs  for  coherent  couplings,  e.g.  a  transverse  field  h  ]T\  aj. 
Our  goal  in  what  follows  is  to  obtain  analytic  expressions  for  such  correlation  functions 
in  full  generality,  and  then  apply  them  to  calculating  a  variety  of  experimentally  relevant 
quantities.  In  order  to  describe  Q  concisely,  it  is  useful  to  define  a  variable  a3  on  each 
site  such  that  a3  —  i  if  there  are  no  occurrences  of  the  operators  aj  in  £,  and  aj  —  0 
otherwise.  Now  we  recognize  that  if  an  operator  <7°=±  occurs  in  Q  one  or  more  times, 
the  operator  dj  (appearing  in  the  time  evolution  operator)  is  forced  to  take  on  a  well 
defined  value  u)(t)  at  all  points  in  time  (see  Fig.  2).  As  a  result,  we  can  rewrite  the 
correlation  function  Q  as 


(  -*  Jc  (lt  *(t)))llfo  +  *i 


(7) 
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Figure  2.  Graphical  representation  of  a  sample  correlation  function  Q  = 
(7c  (0)^J (0 ))  )•  Here  aj  =  =  0,  oii^j.k  =  1,  and  the  time-dependent 

functions  <rj  (t)  and  of.  (t)  are  shown  in  the  bottom  two  panels. 


where 


Pi  =  (i  -  aj)fj  [oj(o)]  fj  kK0*)] »  (8) 

(t  =  0*  marking  the  end  of  the  backwards  trajectory,  Fig.  2),  /  is  the  complex  conjugate 
of  /,  fc  is  a  time  integral  that  runs  along  the  closed-time  path,  and 

£(0  =  \Y1  K1  -  “iK1  -  +  2(!  -  •  (9) 

j,k 


The  first  thing  to  notice  is  that,  since  fc  dt  =  0,  the  final  time-independent  term  in  1C 
vanishes.  Since  this  is  the  only  non  separable  term  in  /C,  the  remaining  time  evolution  is 
straightforward  to  compute.  It  is  helpful  to  define  the  following  parameters  that  depend 
on  the  functions  <Jz-(t) 


Tk  =  J^1 

jyk 


dt  a* ( t ) 


d 


—  2  ^  "  Jjk(l  ai)(l 

j,k 


(10) 

(11) 


in  terms  of  which 

[  dt  lC(t)  =  d  +  ipk°zk- 

Jc  k 

The  time-ordered  correlation  functions  of  Eq.  (7)  can  now  be  compactly  written 

g = e-"  n  m  m  n  (ft + a>) 

i  j 

=  e-n  [ft +*>,)], 


(12) 


(13) 

(14) 


with 


gf(.v)  =  a,  (mWe-*  ±  l/i(-l)|2ei*’)  (15) 

(gj  will  be  useful  momentarily).  The  equivalence  between  Eqs.  (13)  and  (14)  can 
be  understood  by  explicitly  comparing  the  expressions  inside  the  product  for  the  two 
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possible  situations  aj  =  0, 1:  If  oq  =  0,  then  tpj  =  0  and  the  expectation  value  is  unity, 
whereas  when  aj  =  1  we  find  pj  =  0  and  the  expectation  value  gives  g+(pj). 

The  insertion  of  an  operator  crj  (t)  inside  a  correlation  function  Q ,  which  we  denote 
by  writing  Q  — >  Gj,  is  relatively  straightforward.  If  aj  =  0,  then  clearly  the  substitution 
crj  — >  crj(t)  does  the  trick.  If  aj  =  1,  dj  can  be  inserted  by  recognizing  that  the  variable 
Pj  couples  to  dj  as  a  source  term,  and  thus  the  insertion  of  d|(t)  is  equivalent  to  applying 
to  Q.  Both  possibilities  are  captured  by  writing 

Q]  =  ^(1  -  aj)azj{t)  +  aji-^j  G,  (16) 

which,  using  d g+(p)/dp  =  —ig~  (</?),  can  be  simplified  as 

Gj  =  e~ W  [Pj<7j(t)  +  9j  (<Pj)]  n  \?k  +  gk  M]  ■  (17) 

kpj 

Notice  that  if  all  operators  occur  at  the  same  time  t,  e.g.  when  calculating  equal¬ 
time  correlation  functions,  then  '0  —  0  and  pj  =  Jjk{  1  —  ak)aj  (±2 1)  (with  the  ± 
depending  on  whether  d^  is  applied  to  the  spin  on  site  k). 


2.1.  Bloch  vector  dynamics 


It  is  now  straightforward  to  calculate  the  dynamics  of  the  Bloch  vectors 

«i(0  =  i{WW).KW).<^W)}-  (18) 

Because  dj  commutes  with  the  Ising  interaction  the  z  component  of  spin  is  time 
independent,  and  given  by  S?  =  |^(0)  =  ^  cos  0j.  The  transverse  spin  components 
Sj(t)  and  Sj(t)  can  be  obtained  from  the  real  and  imaginary  parts,  respectively,  of 
(d/(t)).  A  straightforward  application  of  Eq.  (13)  gives 


u 


t(t))  =fJ(l)fj(-l)l[g+(2Jjkt) 

1 
2 


=  - e l<t>j  sm 


n  (cos  2  Jjkt  —  i  sin  2Jjkt  cos  9k) . 


(19) 


For  the  special  case  where  all  spins  point  along  9  —  n/2  at  t  —  0  there  is  pure  decay  of 
the  Bloch  vector  without  any  rotation.  In  Fig.  (3)  we  show  the  projection  into  the  xy 
plane  of  the  Bloch  vector  So(t)  (where  j  —  0  labels  the  central  site  of  a  55-site  triangular 
lattice)  for  an  initial  state  in  which  all  spins  point  in  a  single  direction  lying  outside 
of  the  xy  plane  (9  —  7t/4,  0  =  0).  The  Bloch  vector  spirals  inwards:  The  precession 
can  be  understood  as  a  mean-field  effect  [33],  with  the  spin  rotating  due  to  the  average 
magnetization  of  the  other  spins,  while  the  decay  is  due  to  the  development  of  quantum 
correlations.  Note  that,  in  a  finite  system,  the  length  S(t)  =  |<S'(£)|  of  the  total  Bloch 
vector  S(t )  =  JT  Sj(t)  decays  even  at  the  mean- field  level  for  any  (  ^  0.  This  decay, 
however,  is  due  to  the  existence  of  a  spatially  inhomogeneous  mean-field,  and  cannot  be 
associated  with  the  development  of  correlations. 
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Figure  3.  Trajectories  of  the  Bloch  vector  S’o  projected  into  the  xy  plane,  for  all- 
to-all  (left)  and  dipolar  (right)  couplings.  Here  j  —  0  labels  the  central  site  (green 
dot  in  left  pannel)  of  a  55-site  triangular  lattice,  and  all  spins  are  initialized  at 
{0j,4>j}  =  {tt/4,0}.  In  both  cases  we  choose  a  nearest  neighbor  coupling  J.  and  scale 
the  time  by  Jtot  =  -V(rj  —  ro)'“-  This  rescaling  of  time  is  used  so  that  different 
range  interactions  give  rise  to  comparable  precession  rates.  Note  that  at  mean-field 
level  the  trajectory  would  close  on  itself.  The  inward  spiral  indicates  the  growth  of 
quantum  correlations  and  resultant  decay  of  the  spin  length,  and — in  these  rescaled 
time  units — is  more  significant  for  shorter-range  interactions. 


2.2.  Equal  time  correlation  functions 

Spin-spin  correlation  functions  can  be  calculated  just  as  easily  from  Eqs.  (13)  and  (16). 
All  two-point  correlation  functions  can  be  calculated  from  the  four  quantities 


(20) 

c* 

III 

'-•4- 

c-f. 

S 

(21) 

(22) 

=  (<^(0^(0) 

(23) 

and  their  complex  conjugates.  Since  the  Hamiltonian  commutes  with  all  <tJ,  the  first  one 
is  given  trivially  by  Czf  =  gj(0)g^(0)  =  cos  6j  cos  6k ■  The  second  one  can  be  obtained 
from  Eqs.  (16)  and  (19)  as 

c?  =  II  stw) 

l^j,k 

=  sin  Ojgf  (2 Jjkt)  fj  gf  ( 2Jjtt ),  (24) 

and  the  third  and  fourth  are 

Cfk+  =  ^et(4>3+4>k)  sin Oj  sin  6k  ]T[  gf(2Jjit  +  2  Jkit) 

OkJ 

c?k~  =  sin dj  sin  0k  gf{2J]l  -  2  Jklt) 

These  correlation  functions  can  be  used,  for  example,  to  calculate  the  time 
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time  (h/J) 


Figure  4.  Optimal  spin  squeezing  £  (obtained  be  minimizing  £(f)  over  time)  for  a 
variety  of  power-law  couplings  Jij  =  J/\ri  —  r3  \ ^ .  The  different  curves  correspond  to: 
infinite  ranged  ((  =  0,  black  solid  line),  coulombic  ((  =  1,  blue  dashed  line),  dipolar 
(C  =  3,  red  dotted  line),  and  nearest  neighbor  ((  =  oo,  green  dot-dashed  line).  For  this 
calculation  we  take  the  spins  to  all  point  along  the  x-axis  at  t  =  0,  and  use  55  sites  of 
a  triangular  lattice  (the  same  as  shown  in  Fig.  2). 


dependence  of  the  spin  squeezing  parameter 

m  =  (25) 

Here  ASmin(t)  is  defined  to  be  the  minimum  uncertainty  along  a  direction  perpendicular 
to  -S,(f)||.  The  squeezing  parameter  determines  the  phase  sensitivity  in  a  suitably 
performed  Ramsey  experiment,  which  is  enhanced  over  the  standard  quantum  limit 
whenever  £  <  1  [34],  For  (  =  0  (infinite-range  interactions)  the  calculation  was  first 
performed  in  [34],  However,  in  many  experimentally  relevant  situations  the  interactions 
have  some  finite  range  and  the  maximum  achievable  squeezing  is  diminished  (Fig.  4). 

2.3.  Unequal-time  correlation  functions 

It  is  also  possible  to  calculate  correlation  functions  involving  the  application  of  spin 
operators  at  different  times,  which  describe  the  propagation  in  time  of  a  perturbation 
to  the  system.  As  an  example,  we  can  easily  calculate  dynamical  response  functions  of 
the  form 

—  siii6,siri6;s  y+ (2atiJ ^  +  2H2J jk] .  (26) 

These  can  be  combined  to  calculate  dynamical  response  functions  involving  arbitrary 
Pauli  matrices,  some  examples  of  which  are  shown  in  Fig.  5. 

As  we  will  see  in  Sec.  3,  such  dynamical  correlation  functions  allow  us  to  calculate 
the  effect  of  spontaneous  relaxation  and  excitation  on  the  dynamics  of  the  system. 

||  If  we  choose  our  x-axis  to  be  in  the  direction  of  S(t),  and  define  =  \  JT  (cos(^)i7|  +  sin (b)crj), 
then  A5min(t)  is  obtained  by  minimizing  (( S —  (S^)2)1^2  over  if). 
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(b) 


0  7 r 

t  (units  of  h/J) 


For  instance, 
relaxation  of 
given  by 


Figure  5.  Connected  dynamical  correlation  function  (t\ ,  t2)  =  ,  t-2)  — 

for  a  100  site  ID  chain  with  £  =  1  and  nearest-neighbor  coupling  J. 
(aVwh  plot  S#+r(0,i)  for  i  =  50  and  r  =  {2,  3,4,5}  (from  top  to  bottom).  In  (b) 
p/ot  $\Vi+ 1(t,  t  +  St)  (again  for  i  =  50)  as  a  function  of  t  and  St.  In  both  plots  the 
ne  consists  of  all  spins  pointing  in  the  +x  direction  (6  =  7t/2  and  cf>  =  0). 


injhratlj^olarized  along  the  x  axis,  the  effect  of  the  sudden 
a  spin  on  site  k  at  time  ts  on  the  time  dependence  of  er^  (for  j  ^  k )  is 


(i s)>  =Re 


^2 ijjk(t  2 ts) 


]  J  cos(2  Jjit) 
=  cos(2J,tf[l  -  2 tjt})  {8j(t))  . 


(27) 


Note  that  this  is  the  same  time  evolution  we  would  obtain  if  no  spontaneous  relaxation 
had  occurred,  the  fc’th  spin  were  simply  absent,  and  the  j'th  spin  were  coupled  to  a 
longitudinal  magnetic  held  of  strength  2Jjk(2ts/t  —  1).  The  reason  for  this  behavior  is 
straightforward  when  considering  the  time  evolution  in  the  Schrodinger  picture.  The 
application  of  to  the  wave  function  at  time  ts  not  only  forces  the  fc’th  spin  to  point 
down  between  times  ts  and  t,  but  also  destroys  the  piece  of  the  initial  wave  function 
having  weight  into  states  with  <Jk  =  —1.  Hence  it  is  as  if  the  fc’th  spin  pointed  up  for  a 
time  tUp  =  t8,  and  down  for  a  time  tdown  =  t  —  ts,  thus  contributing  an  inhomogeneous 
longitudinal  magnetic  held  of  strength  2 Jjk(tup  —  tdown)/^  =  2Jjfc(2ts/f  —  1)  (the  factor 
of  2  arises  because  the  ,/jfc  couple  to  Pauli  matrices  rather  than  the  spin-f  matrices). 


3.  Inclusion  of  dissipation 

In  the  previous  sections  we  have  treated  our  Ising  spins  as  a  closed  system.  That  is, 
we  neglected  any  coupling  that  might  exist  between  our  system  and  the  outside  world, 
and  thus  initially  pure  states  remained  pure  throughout  the  dynamics.  In  any  physical 
realization  of  the  Ising  model,  this  is  clearly  an  idealization;  decoherence  occurs  and  often 
must  be  accounted  for.  For  example,  Rydberg  atoms  suffer  from  spontaneous  emission 
[53],  while  ions  couple  strongly  to  fluctuating  classical  (electric  and  magnetic)  helds 
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and  can  decohere  through  off-resonant  light  scattering  from  the  spin-dependent  optical 
dipole  forces  used  to  engineer  the  Ising  interactions  (this  off-resonant  light  scattering  can 
produce  spontaneous  excitation/relaxation  and  dephasing)  [54],  One  way  to  envision 
dynamics  in  an  open  system  is  by  considering  the  probabilistic  occurrence  of  sudden 
perturbations  of  the  system — quantum  jumps — due  to  the  system-environment  coupling 

[55] .  As  suggested  in  Sec.  2.3,  and  as  will  be  explained  in  detail  below,  the  strategy 
we  have  developed  for  computing  unequal-time  correlation  functions  is  well  suited  to 
describing  such  effects. 

3.1.  Description  of  the  problem 

Given  a  density  matrix  g  describing  a  system  coupled  to  a  reservoir,  it  is  always  possible 
to  express  the  expectation  value  of  a  system  operator  A  in  terms  of  the  system  reduced 
density  matrix  p  =  trR  [p]  as  (A)  =  trgfpM].  Here  trR(s)  denotes  a  trace  over  the  reservoir 
(system)  degrees  of  freedom — we  will  drop  these  subscripts  from  now  on,  since  all  future 
instances  of  tr  refer  to  a  trace  over  system  degrees  of  freedom  only.  In  this  language, 
the  effect  of  a  finite  system-reservoir  coupling  is  that  an  initially  pure  system  density 
matrix  p(0)  =  |\Er(0))(\Er(0)|  will  evolve  into  a  mixed  state  (reflecting  entanglement 
between  the  system  and  reservoir  degrees  of  freedom).  When  the  system-environment 
coupling  is  weak  (i.e.  small  compared  to  the  inverse  of  relevant  system  time  scales)  and 
the  reservoir  correlation  time  is  small,  the  Born-Markov  approximation  is  justified  and 
the  reduced  system  density  matrix  obeys  a  Markovian  master  equation  of  Lindblad  form 

[56] .  We  choose  a  very  general  master  equation  appropriate  for  describing  the  various 
types  of  decoherence  relevant  to  trapped  ions  [54],  Rydberg  atoms  [53],  and  condensed 
matter  systems  such  as  quantum  dots  [5]  and  nitrogen  vacancy  centers  [6]: 


p  =  -iJf?(p)  -  «^ud(p)  -  ^du(p)  -  2S?ei (p),  (28) 

where 

JT(p)  =  [ %p]  (29) 

^ud(p)  =  NdJ p  +  Pata7  ~  2a7 Pat )  (30) 

3 

^du(p)  =  %  p  +  p*7 ~  PaJ )  (31) 

3 

-^ei(p)  =  yE(2/,“  2&jP&j)  ■  (32) 

3 


The  first  term  involving  a  commutator  describes  coherent  evolution  due  to  the  Ising 
interaction,  and  the  various  terms  having  subscripts  “ud”,  “du”,  and  “el”  correspond 
respectively  to  spontaneous  relaxation,  spontaneous  excitation,  and  dephasing^f. 


The  “ud”  and  “du”  subscripts  remind  us  that  the  corresponding  decoherence  processes  change  a  spin 
state  from  “up”  to  “down”  (or  vice  versa)  along  the  z  axis,  while  the  “el”  subscript  reminds  us  that 
dephasing  is  “elastic”  in  the  sense  that  it  does  not  change  the  spin  projection  along  the  2  axis. 
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Equation  (28)  has  the  formal  solution  p(t)  =  ^(t)p( 0),  with 

(t)  =  exp  [- 1  (iAf  +  =Sfud  +  Xiu  +  =£?ei)]  •  (33) 


The  exponential  of  super-operators  is  meant  to  be  understood  via  its  series  expansion, 
in  which  the  multiplication  of  two  Lindblad  super-operators  implies  composition 
(NT\  xNT2)  (p)  =  (J2?2  (p))-  Our  goal  in  what  follows  is  to  compute  the  time  dependence 

of  an  arbitrary  operator  A  at  time  t,  given  in  the  Schrodinger  picture  by 


A{t)  =  tr  Ap{t) 


(34) 


3.2.  Dephasing  (T2  decoherence) 


An  immediate  simplification  follows  from  the  observation  that 


—  [Jzfei,  Jfud]  —  0. 


(35) 


That  the  last  two  commutators  vanish  is  less  obvious  than  the  first,  but  physically  it 
has  a  very  clear  meaning:  Spontaneous  relaxation/excitation  on  a  site  j  causes  the  j' th 
spin  to  have  a  well  defined  value  of  crj,  and  thus  to  be  unentangled  with  the  rest  of  the 
system.  Since  the  dephasing  jump  operator  cfj  changes  the  relative  phase  between  the 
states  \<jj  =  ±1),  whether  spontaneous  relaxation/excitation  occurs  before  or  after  a 
dephasing  event  only  affects  the  sign  of  the  overall  wave  function,  which  is  irrelevant. 
As  a  result,  we  can  write 


%  (p)  =  e~ +^ud+-2du) 


(36) 


and  the  time  dependence  of  an  arbitrary  observable  A(t)  =  Tr 


p{t)A 


can  be  written4" 


A(t)  =  Tr  e  h*^+^ud-K2du)p(Q)e 


(37) 


Note  that  application  of  a  superoperator  does  not  commute  with  operator  products 
(Jf(01)02  ^  Jf(0i02)),  and  we  use  the  convention  that  a  superoperator  should  act  on 
the  operator  immediately  to  its  right.  The  effect  of  the  time  evolution  due  to  JCe\  can 
be  understood  by  considering  its  effect  on  the  Pauli  operators: 


e-^  a]’y  =  e-r^t/2axfy, 


=  d‘. 


(38) 


In  light  of  equations  (37)  and  (38),  we  are  free  to  ignore  the  dephasing  terms  in  the 
master  equation  at  the  expense  of  attaching  a  factor  of  e~Telt /2  to  every  operator  ax-'y 
occurring  inside  an  expectation  value: 


Tr 


p(t)A(ax,& 


->•  Tr 


p(t)A(e 


-relt/2~x  e-relt/2^y 


ai  > e 


j  ■ 


(39) 


where  p(t)  on  the  right-hand  side  evolves  under  the  master  equation  without  the 
dephasing  term. 


Note  that  we  apply  the  operator  e  *-s?el  to  A  rather  than  p(0).  This  is  justified  by  the  identity 


tr 


^el(Ol)C?2 


=  tr 


0^61(02) 


,  true  for  arbitrary  operators  0\  and  02l  which  holds  because  the 


jump  operators  a J  are  hermitian. 
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3.3.  Spontaneous  relaxation  and  excitation  (T\  decoherence) 

Because  the  effects  of  dephasing  are  fully  included  by  Eq.  (39),  the  remaining  problem 
is  to  compute  the  time  dependence  of  operators  whose  expectation  values  are  taken  in 
a  density  matrix  evolving  simultaneously  under  Jif,  J2?ud,  and  2%u: 

A(t)  =  Tr  e-t(i^+^ud+^du) (40) 

Formally  the  challenge  of  including  the  effects  of  spontaneous  relaxation  and  excitation 
is  related  to  the  nontrivial  commutation  relation 

[jr,^ud(du)]  T^o.  (4i) 

Physically,  the  obstacle  is  that  spontaneous  relaxation  and  excitation  change  the  value 
of  <t~  for  the  spin  which  they  affect,  and  this  change  feeds  back  on  the  system  through 
the  Ising  couplings.  From  the  results  on  coherent  dynamics  presented  in  Sec.  2,  we 
know  that  time  evolution  under  alone  is  tractable.  This  suggests  that  we  attempt 
to  solve  Eq.  (40)  by  rewriting 

iAif  +  d  +  =  idrffe  ff  —  (42) 

where  M  contains  all  and  only  terms  that  do  not  commute  with  AC ,  and  then  doing 
perturbation  theory  in  M.  The  above  separation  is  accomplished  by  defining  an  effective 
(non-Hermitian)  Hamiltonian  and  its  corresponding  superoperator 

HeS  =  iJV1^  and  ACeS(p)  =  W-cEp  ~  pU]eS,  (43) 

j 

and  the  recycling  term 

& (P)  =  rud  +  rdu  P°J  ’  (44) 

j  j 

in  terms  of  which  the  time  evolution  operator  is  W(t)  =  In  Eq.  (43)  we 

have  defined  7  =  |(Tud  —  Tdu)  and  Tr  =  Tud  +  Tdu-  Defining  %)(t)  =  e~ztJl)feff ,  we  can 
now  expand  the  time  evolution  operator  as  a  power  series  in  &  in  order  to  obtain  the 
time-dependent  expectation  value  A(t): 

Ait)  = 

rt  rt.2  P 

y,  /  dtn  .. .  dt  1  tr  A%(t  -  tn)&%(tn  -  tn- 1) . . .  %(t2  -  ti)&%(ti)p(0)  .  (45) 

n  J° 

This  expansion  is  the  underlying  object  being  evaluated  when  Monte  Carlo  wave  function 
methods  [55]  (quantum  trajectories)  are  used  to  approximate  the  density  matrix.  In 
Appendix  A,  we  show  in  detail  how  the  series  in  Eq.  (45)  leads  to  an  expression  for 
A(t)  in  terms  of  the  closed-time  path  ordered  correlation  functions  obtained  in  Sec.  2, 
and  the  summation  of  that  series  is  carried  out  in  Appendix  B.  Here  we  will  simply 
summarize  the  calculation,  and  explain  in  physical  terms  the  essential  structure  of  the 
Hamiltonian  and  decoherence  that  allows  the  result  to  be  cast  in  closed  form. 

We  begin  by  noting  that  when  writing  A{t)  as  a  sum  over  closed-time  path  ordered 
correlation  functions,  each  operator  inserted  along  the  forward  leg  of  the  time-contour  is 
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accompanied  by  its  hermitian  conjugate  appearing  at  the  same  time  along  the  backward 
part  of  the  time  contour.  If  we  had  explicitly  included  an  environment  and  attempted 
to  trace  over  it  (rather  than  starting  with  a  Markovian  master  equation),  this  feature 
of  the  problem  would  emerge  as  a  direct  consequence  of  the  Markov  approximation. 
As  a  result,  it  suffices  to  describe  any  term  in  the  series  expansion  by  specifying  the 
occurrence  of  operators  on  the  forward  time  contour.  To  facilitate  this  description,  we 
introduce  notation  describing  the  occurrence  of  operators  belonging  to  a  particular  site 
j  (see  Fig.  6  for  a  summary).  We  take  TZf  to  be  the  number  of  times  the  operator  of 
occurs  along  the  forward  time  path,  ,t^  }  to  be  the  set  of  times  at  which  jump 

operators  are  applied  to  site  j,  TZj  =  7 Zf  +  7 Zj,  Kj  =  ±1  depending  on  whether  the 
operator  at  the  latest  time  along  the  forward  path  is  of,  and 

Tj  =  (1  -otj)  [  (46) 

Jo 

We  will  also  use  bold  symbols  7Z,  k,  and  r  to  specify  the  complete  set  of  these  variables 
on  all  lattice  sites.  Note  that  specifying  7 Zj  and  Kj  determines  both  IZj  and  7 Zj , 
since  two  consecutive  (in  the  closed-time-path-ordering)  applications  of  an  operator  of 
gives  zero.  Therefore,  we  will  only  include  the  7 Zj  and  Kj  as  explicit  arguments  in  the 
correlation  functions  below.  These  variables  are  sufficient  to  determine  the  value  of  any 
term  in  the  series  expansion  of  A(t),  so  we  do  not  need  to  keep  track  of  the  individual 
times  at  which  each  jump  operator  is  applied. 
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Figure  6.  A  few  examples  of  how  spin-raising  and  spin-lowering  operators  belonging 
to  the  j ’tli  site  may  occur  along  the  forward  time  evolution,  and  the  notation  used  to 
characterize  these  occurrences. 

All  nonzero  terms  in  Eq.  (45)  are  captured  by  summing  over  the  IZj  and  Kj,  and 
integrating  over  the  times  t[ , . . . ,  ,  denoted 
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If  A  can  be  written  as  a  product  of  operators  ab>  ( b3  =  ±)  on  sites  j  contained  in  a  set 
77,  then  A{t)  can  be  compactly  expressed  (see  Appendix  A)  as 

A(t)  =  j  ^  r)g(s,  n,  r).  (48) 

The  prefactor 

T(U,  «,  r)  =  ll  VjiUj,  Kj,  Tj)  =  J]  (e-W  (Tud)^"  (Tdu)^+  e"2^')  ,  (49) 

j  j 

is  closely  related  to  the  probability  that  a  series  of  jumps  described  by  the  variables  n, 
k,  and  r  has  occurred.  Defining  the  symbol  s  as  the  vector  of  quantities  s3  =  (pj/t  —  2i'y, 
the  correlation  functions 

g(s,n,T)  =  Yl^rEr  Tj)  =  (er*WTi)  [Pj  +  gfisjt)])  ,  (50) 

j  j 

are  of  the  general  form  presented  in  Sec.  2.  Careful  bookkeeping  reveals  that  the 
variables  ip  and  $  defined  in  Sec.  2  are  given  by 


Tj 

^  ^  bkJjki 

(51) 

ker) 

d(r) 

= 

j 

(52) 

0j(Tj) 

= /  ° 

[  2  Tj'Tllk&r1^kJjk 

if  j  e  77; 

if  j  £  77. 

(53) 

Note  that  the  7Z3  dependence  in  Qj  is  hidden  in  the  implicit  dependence  of  pj  and  on 
o/,3  (which  for  j  ^  rj  satisfies  a3  =  Sn:i,o)-  When  evaluating  the  sums  and  integrals  in  Eq. 
(48),  one  must  keep  in  mind  that,  whenever  j  G  77,  terms  with  1Z3  ^  0  vanish  because 
they  contain  the  consecutive  application  of  either  cr+  or  a~ .  Physically,  the  vanishing  of 
such  terms  reflects  the  lack  of  coherence  for  any  spin  that  has  undergone  even  a  single 
spontaneous  spin  flip. 

The  factorization  of  V  into  functions  Vj  of  local  site  variables  {TLV  k3,  t3}  is  a  direct 
consequence  of  the  single  particle  nature  of  the  anti-Hermitian  part  of  TLes.  Physically, 
this  factorization  occurs  because  the  dissipation  we  are  considering  is  uncorrelated  from 
site  to  site  (in  contrast  to  the  collective  relaxation  processes  that  arise,  e.g.,  in  the 
context  of  Dicke  superradience  [57]).  The  factorization  of  the  correlation  function 
Q  into  functions  Qj  of  local  site  variables  {7Z3,t3}  is  a  more  surprising  result,  and 
depends  crucially  on  the  occurrence  of  jump  operators  at  the  same  times  on  the  forward 
and  backward  time  evolution  (Appendix  A).  The  7  appearing  in  the  argument  of 
dt  =  Tj  ~~  2 ijt)  affects  the  value  of  any  term  in  Eq.  (48)  in  which  no  jump  operators 
are  applied  to  site  j  (such  that  aj  =  1  and  therefore  (s3t)  7^  0).  It  arises  from  the 
term  — iqd2  in  TLes  [Eq.  (43)],  and  decreases  the  expectation  value  of  <jJ  for  7  >  0  (when 
spontaneous  relaxation  outweighs  spontaneous  excitation).  This  effect  is  often  referred 
to  as  null  measurement  state  reduction :  Gaining  knowledge  that  the  spin  on  site  j  has 
not  spontaneously  flipped  affects  the  expected  value  when  measuring  <jJ. 
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Because  Q  and  V  both  factorize,  we  need  only  to  evaluate  the  quantities 

(sj >  * )  =  f  Vi ’  Ko ’  rt ) $3  (si >  nj Nj)  (54) 


in  terms  of  which 


'A{t)  =  Y[$j{sj,t).  (55) 

3 

The  explicit  dependence  on  Sj  is  included  to  remind  the  reader  that,  after  the  sums  and 
integrals  (over  7 Zj,  kv  and  t\  . . .  tJni)  have  been  carried  out,  Sj  is  the  only  site-dependent 
quantity  on  which  $j(sj,t)  depends.  We  evaluate  these  sums  and  integrals  in  Appendix 
B,  obtaining 


t ) 


'  e-vvt/2 


Pj 


=-rr/2 


cos  ft Ns‘j  —  r  j  +  +  ipj  cos  9j)  t  sine  ft  Cs‘j  —  r 


if  j  e  r/; 

d  (56) 
if  J  f  V, 


where  r  =  Tudrdu. 

If  A  also  contains  an  operator  erf  (t)  (with  /  ^  r/),  we  denote  its  expectation  value 
by  Af(t).  The  insertion  of  af  must  be  dealt  with  at  the  point  of  Eq.  (54).  Keeping  in 
mind  the  discussion  surrounding  Eq.  (16),  and  remembering  that  Si  =  ipi/t  —  2iy,  we 
must  replace  &i(si,t)  with 

Vi{si,t)  =  y ^(1  -  ai)Ki  +  ail--^j  Gi(shHhTi)  (57) 
Therefore  we  have 


Al(t)  =  ^i(sht)Y[^j(sj,t), 
and  in  Appendix  B  we  find 


'h i(sht )  =  e 


_  p-rr/2 


cos  (t  yjs'j  -  r  )  +  (isi  +  27  -  y  cos  6i  \  t  sine  (t  y4 ?  -  r  j 


(58) 


(59) 


3-4-  A  simple  application:  under-damped  to  over-damped  transitions 


These  equations  reveal  that  correlation  functions  will  generally  undergo  a  qualitative 
transition  in  dynamics — from  over-damped  to  oscillatory — whenever  the  condition  sf  = 
r  is  satisfied.  This  behavior  is  the  most  clearly  manifest  when  the  couplings  Jt]  have 
a  simple  structure,  such  as  nearest  neighbor  or  all-to-all.  For  instance,  for  nearest- 
neighbor  coupling  in  ID,  assuming  {Tei,  Tud,  Tdu}  =  {0,T,T},  and  choosing  the  initial 
state  to  point  along  the  x— axis,  we  find  (ignoring  boundary  effects) 


Sx(t )  =  —e~3rt 
y  J  2 


cos 


tV  4  j2  -  r2 )  +  rt  sine  ( tV4 J2  -  r2 


(60) 


which  becomes  critically  damped  at  Tc  =  2  J  (see  Fig.  7).  It  is  interesting  to  note  that 
this  solution  is  similar  in  structure  to  the  damping  of  a  classical  harmonic  oscillator 
or  a  coherently  driven  two-level  system  (c.f.  the  weak-coupling  limit  of  the  spin  boson 
problem  [58,  14]).  It  is  important  to  contrast  this  behavior  with  that  of  a  single  spin 
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Figure  7.  Dynamics  of  the  transverse  spin  length  Sx  in  a  ID  nearest  neighbor  Ising 
model,  with  {rei,  ru(j,  Tdu}  =  {0,r,r},  for  several  different  values  of  T  between  0 
(undamped)  and  Fc  (critically  damped):  T  =  0  (black  solid  line),  T  =  Tc/8  (blue 
dashed  line),  F  =  Tc/4  (red  dotted  line),  and  T  =  Fc  (green  dot-daslred  line). 


coupled  to  a  Markovian  bath,  where  the  decoherence  we  consider  only  causes  a  damped- 
to-oscillatory  transition  in  the  presence  of  a  transverse  magnetic  field;  the  Hamiltonian 
dynamics  must  be  able  to  restore  coherence  in  the  basis  for  which  the  environment 
induces  a  measurement.  In  the  present  case,  there  is  no  transverse  field,  and  it  is  not  a 
priori  obvious  that  such  behavior  should  emerge.  In  fact,  a  simple  mean-field  estimate 
of  the  dynamics  fails  to  capture  the  oscillatory-to-damped  transition.  Using  a  site- 
factorized  ansatz  for  the  density  matrix,  p  =  0  pj,  it  is  straightforward  to  see  that 
[33] 

Smf(C  =  ye~rtcos(4Jcos(#)i).  (61) 

Thus  mean-field  theory,  which  assumes  an  unentangled  density  matrix,  always  predicts 
under-damped  dynamics;  the  transition  to  over-damped  behavior  captured  by  the  exact 
solution  depends  crucially  on  the  competition  between  decoherence  and  entanglement. 

3.5.  Qualitative  insights  into  decoherence  in  interacting  many-body  systems 

In  addition  to  providing  an  efficient  way  to  compute  arbitrary  observables  for  an  open 
many-body  system,  the  above  calculation  provides  significant  insight  into  the  interplay 
beween  decoherence  and  interactions.  To  keep  the  notation  as  simple  as  possible,  we  will 
focus  on  the  case  of  nearest-neighbor  coupling  on  a  lattice  with  coordination  number 
z,  and  calculate  the  time  dependence  of  the  transverse  spin  length  of  a  single  spin 
sx(z,t )  =  ( <jx )  (for  an  infinite  system  it  is  independent  of  j)  for  a  state  in  which  all 
spins  are  initially  polarized  along  the  x-axis.  In  the  absence  of  decoherence  but  in  the 
presence  of  a  longitudinal  magnetic  field  of  strength  h ,  application  of  results  in  Sec.  2 
gives 

sxh(z,t)  =  Ue  [e-2i“cos‘(2Ji)]  . 


(62) 
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In  the  absence  of  a  longitudinal  field  but  including  equal  rates  of  spontaneous 
relaxation/excitation  (Tudjdu  =  T,  7  =  0),  we  find  instead  [from  Eqs.  (55)  and  (56)] 


sx(z,t )  =  e  b+1)rT-i/2  cos  (t\/AJ'2  —  r2)  +  R sine  (fV4J2  —  T2 


(63) 


However,  it  is  instructive  to  temporarily  hold  off  evaluating  the  sums  and  integrals 
implicit  in  the  <hj,  and  work  directly  with  Eqs.  (54)  and  (55): 

NN  D 

sx(z,t )  =  $j{sj,t)  JJ$fc(sfc,f)  JJ$i(sj,t).  (64) 

k  l 

In  Eq.  64  we  have  divided  the  product  over  lattice  sites  into  three  parts:  the  j’th  site, 
the  nearest  neighbor  sites  (product  labeled  by  NN)  and  the  rest  of  the  lattice  (product 
labeled  by  D,  reminding  us  that  these  sites  are  Disconnected  from  site  j).  First,  we 
observe  that  <3 ?j(sj,t)  =  ^e~Flt^2.  Next  we  evaluate  nf  ^(sz  —  0,  t)  —  1  [the  si  =  0 
because  these  sites  are  disconnected  from  the  j’th  site,  and  hence  Jji  =  0,  see  Eq. 
(51)],  which  follows  from  a  sum  rule  77(77;,  Kh  Ti)Gi(0,  77;,  t{)  =  1  (derivable  from 
tr  [p(t)\  =  1).  It  remains  only  to  evaluate 


NN  „  /NN 

E/  (n^, 


KjTj)Gj(sj 


2  j,nj,Tj) 


(65) 


where  we’ve  adopted  an  abbreviated  notation  ]T)NN/'  =  n^N  f°r  the  sums  and 
integrals  over  the  nearest-neighbor  sites.  Utilizing* 

NN 

nSjPJ.Wf.Tj)  =  2-Ke-2"“  cos(2  Jty~K  =  —s&z-K't),  (66) 

j 

where  7Z  =  7Z3,  t  =  ,  and  h  =  Jr/t ,  and  noting  that  s^(z  —  7 Z,t)  only 

depends  on  r  and  7 Z  (and  not  any  other  combinations  of  the  local  site  variables),  we 
can  then  write 


Z 


\z,t)  =  e-r-t/2  ^ 

TZ=0 


dh  P(n ,  h)sxh(z-K,t)- 


(67) 


Here 

NN  /NN  ,  \ 

P(W,  ft')  =  A  Y.J  (n  CN0H1  j  S(T  _  T')5R,R,  (68) 

is  obtained  by  carrying  out  the  sums  and  integrals  while  holding  77  and  r  fixed  (the 
factor  of  t/J  arises  when  changing  variables  from  r  to  h  in  the  remaining  integral). 

Equation  (67)  has  a  very  suggestive  form.  The  factor  of  e~Tlt /2  out  front  is  the 
single-particle  contribution  of  Tf  decoherence,  and  would  be  present  in  the  absence 
of  interactions;  it  reflects  the  probability  that  the  j’th  spin  has  not  spontaneously 
flipped  before  the  time  t  (a  prerequisite  for  having  any  coherence  along  x).  As  the 


*  Note  that  in  the  final  equality  of  Eq.  (66),  we  have  assumed  that  the  left-hand  side  is  real.  While 
this  is  not  strictly  true,  the  imaginary  part  of  the  left-hand  side  will  vanish  after  the  integral  over  h  is 
carried  out  below,  so  we  make  no  mistake  by  ignoring  it. 
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Figure  8.  The  interplay  between  Ising  interactions  and  spontaneous  spin  flips  induces 
a  source  of  decoherence  beyond  the  direct  action  of  the  spin  flips  themselves.  In  (a),  the 
j’th  (central)  spin  evolves  due  to  the  Ising  couplings  with  its  neighbors.  Spontaneous 
relaxation/excitation  by  its  neighbors  couples  back  (via  the  Ising  interactions)  as  a 
temporally  fluctuating  longitudinal  field  (b),  inducing  a  non-local  dephasing  process. 


z  nearest  neighbors  evolve  in  time,  they  can  undergo  spin  flips  which  cause  them  to 
fluctuate  in  time  between  pointing  along  +z  and  —2  [Fig.  8(a)],  Even  once  flipped, 
they  influence  (via  the  Ising  couplings)  the  j’th  spin  in  a  manner  formally  equivalent 
to  a  longitudinal  magnetic  field  of  strength  h  =  J(j/t )  [Fig.  8(b)],  The  quantity 
P(  1Z,  h )  describes  the  probability  that  7Z  nearest  neighbors  have  spontaneously  flipped, 
collectively  contributing  an  effective  magnetic  field  of  strength  h  to  the  time  evolution  of 
the  j’th  spin.  The  sum  over  1Z  and  integral  over  h  then  average  the  resulting  dynamics 
for  the  j’th  spin,  s^(z  —  7 Z,t),  over  the  possible  behaviors  of  its  neighbors.  For  a  given 
TZ,  the  integral  over  h  reduces  the  Bloch  vector  length  by  an  amount  depending  on 
the  width  in  h  of  the  distribution  P(1Z,  h).  Physically,  this  integral  captures  the  phase 
diffusion  of  the  j’th  spin  due  to  the  stochastic  temporal  fluctuation  of  its  immediate 
environment  (neighboring  spins,  as  shown  in  Fig.  8).  Thus  we  see  very  clearly  that  the 
interplay  between  spontaneous  spin  flips  and  coherent  interactions  leads  to  an  emergent 
source  of  dephasing:  Flipping  spins  act  as  fluctuating  magnetic  fields  (mediated  by 
the  Ising  interactions)  on  other  spins,  even  if  these  latter  spins  have  not  been  directly 
affected  by  decoherence. 

4.  One-axis  twisting  in  an  open  system 

The  expressions  in  Eqs.  (56)  and  (59)  furnish  a  complete  description  of  correlation 
functions,  and  in  special  cases  afford  descriptions  of  common  experimental  observables 
and  entanglement  witnesses.  As  a  concrete  example,  we  will  use  these  expressions  to 
study  the  development  and  loss  of  entanglement  in  an  open-system  version  of  the  one- 
axis  twisting  (OAT)  model.  It  is  important  to  keep  in  mind,  however,  that  most  of 
the  following  results  can  be  generalized  to  take  into  account  arbitrary  Ising  couplings. 
Defining  the  collective  spin  operator  Sz  =  ^  ]Tb  <rj ,  the  OAT  Hamiltonian  is  given  by 

U  =  2  J(SZ)2,  (69) 
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(a) 


(b) 


Figure  9.  Spin  squeezing  (dB)  of  N  =  103  spins  with  no  decoherence  (black  solid 
line),  pure  dephasing  ({rei,  rud,  rdu}  =  {T,  0,  0},  blue  dashed  line),  and  equal  amounts 
of  spontaneous  relaxation/excitation  ({Tei,  rud,  Tdu}  =  {0,r/2,F/2},  red  dotted  line). 
In  (a)  we  show  the  optimal  squeezing  (optimized  over  angle  ip)  as  a  function  of  time. 
In  (b)  we  plot  the  normalized  variance,  evaluated  at  time  ta,  as  a  function  of  the 
squeezing  angle  ip.  Note  that  at  the  angle  of  maximal  squeezing  (ip  k  0),  dephasing 
is  much  less  detrimental  than  spontaneous  relaxation/excitation,  whereas  both  have 
a  similar  effect  at  the  angle  of  maximal  antisqueezing  (ip  ss  7r/2).  In  all  plots,  we 
choose  r  =  l/ts,  with  ts  =  IbV_2/,3(31/6/2 J)  being  the  time  of  optimal  squeezing  in 
the  absence  of  decoherence. 


which  is  the  (  =  0  limit  of  our  more  general  Ising  model  [Eq.  (1)].  For  an  initial 
state  polarized  along  the  x-axis  ( 6  =  7t/2,  (p  =  0),  it  is  well  known  [34]  that  the  OAT 
Hamiltonian  generates  spin  squeezed  states  at  short  times  [the  squeezing  is  optimal 
at  ts  =  frAf~2,/3(31/6 /2J)],  a  fact  that  has  been  exploited  in  a  number  of  beautiful 
experiments  [35,  36].  In  principle  (i.e.  in  the  absence  of  any  decoherence  or  other 
imperfections),  these  spin  squeezed  states  allow  for  precision  metrology  with  a  phase 
sensitivity  that  scales  as  A/"-5//6,  thus  beating  the  J\f _1/2  scaling  of  the  standard  quantum 
limit.  At  time  t*  =  hn/4:J,  the  OAT  Hamiltonian  gives  rise  to  a  GHZ  (or  Schrodinger 
cat)  state  [59],  which  in  principle  affords  Heisenberg  limited  (~  A^_1)  sensitivity  in  phase 
estimation.  In  the  following  subsections,  we  use  the  results  of  Sec.  3  to  extend  calculate 
spin  squeezing  and  characterize  the  metrological  utility  of  (and  GHZ-type  entanglement 
of)  the  state  at  t*,  which  would  be  the  GHZ  state  in  the  absence  of  decoherence. 

4-T  Spin  Squeezing 

Given  the  results  in  Sec.  3,  analytic  calculation  of  the  squeezing  parameter  in  the 
presence  of  arbitrary  decoherence  rates  Tel,  Tud,  and  Tdu  is  now  straightforward.  As 
can  be  seen  in  Fig.  9(a),  the  effect  of  decoherence  is  much  less  severe  than  that 
of  Ti  decoherence.  One  reason  for  this  behavior  is  that  the  minimum  variance  AS,2lin 
occurs  at  an  angle  ip  (in  the  yz  plane)  only  slightly  deviating  from  the  z  axis.  Therefore 
dephasing,  which  can  be  thought  of  as  random  rotations  of  the  individual  spins  around 
the  z-axis,  does  not  introduce  much  noise  in  the  squeezed  quadrature  (see  Fig.  9b).  To 
the  contrary,  spontaneous  relaxation/excitation  processes  introduce  noise  directly  into 
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Figure  10.  Coherence  of  a  GHZ  state  created  in  the  presence  of  various  types 
of  decoherence:  {Fei,  rud,  rdu}  =  {r,  0,0}  (dephasing,  dashed  blue  line)  and 
{rel,rud,rdu}  =  {0,r,0}  (spontaneous  relaxation,  dotted  red  line).  The  region  above 
the  solid  black  line  is  guaranteed  to  have  TV"-particle  GHZ  type  entanglement. 


the  squeezed  quadrature. 


4-2.  Macroscopic  superposition  states 


In  the  absence  of  decoherence,  the  OAT  Hamiltonian  is  known  to  give  rise  to  A/"-spin 
GHZ  states)} 


|GHZ)  =  1  ^  +^+  I  ^  ,  (70) 

at  a  time  t*  =  Htt/AJ,  where  |  'f}x)  (|  Ij-^))  denotes  the  state  where  all  spins  point  along 
the  positive  (negative)  x-axis  [59].  These  entangled  states  afford  Heisenberg- limited 
phase  sensitivity  [60],  and  are  a  resource  for  certain  types  of  fault-tolerant  quantum 
computation  ([37]  and  references  therein).  However,  they  are  also  a  canonical  example 
of  a  fragile  quantum  state,  and  their  usefulness  is  easily  destroyed  by  decoherence  [61]. 
The  effect  of  dephasing  on  the  production  of  GHZ  state  via  one-axis  twisting  is  well 
understood  [61,  62],  With  the  results  of  Sec.  3,  however,  we  can  easily  calculate 
the  effects  of  dephasing  and  spontaneous  relaxation/excitation  on  the  production  of 
a  GHZ  state  by  one-axis  twisting  in  a  unified  way.  In  this  section  we  explicitly  compare 
the  effects  of  dephasing  to  the  those  of  pure  spontaneous  relaxation  ({Tel,  Tud,  Tdu}  = 
{0,  T,  0}). 


We  first  characterize  the  GHZ  state  by  its  phase  coherence,  obtained  from  the 
expectation  value  C  =  tr  p{t)C  of  the  operator 


c  =  |  i’gX'IJ'z  I 
=^n  (ai+aj-aj)- 


(U) 


0  Strictly  speaking  the  form  given  only  applies  when  the  particle  number  TV"  is  even.  For  odd  TV"  the 
GHZ  state  created  looks  similar,  but  is  composed  of  states  polarized  along  the  ±y  direction. 
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The  quantity  C  characterizes  the  extent  to  which  the  superposition  between  the 
macroscopically  distinct  states  |  'f|\r)  and  |  fj.x)  is  quantum  mechanical  (rather  than 
a  classical  mixture).  Formally,  C  serves  as  a  witness  to  A/"-particle  entanglement  of 
the  GHZ  type,  with  entanglement  guarantied  whenever  \C\  >  1/4  is  satisfied  [63,  37]. 
Application  of  the  results  in  Sec.  3  yields 

C«  =  w  £  £  O  (?)  (2Jt2n  -  "*1  -  »T.  C~m  ,  (72) 

and  in  the  absence  of  decoherence  one  finds  \C(t*)\  =  1/2.  As  can  be  seen  in  Fig.  10,  the 
effect  of  spontaneous  relaxation  on  the  coherence  C  is  comparable  to  (but  worse  than) 
the  effect  of  dephasing.  Both  types  of  decoherence  cause  a  loss  of  phase  coherence  when 
T  ~  J/Af,  with  the  factor  of  Af  responsible  for  the  fragility  of  a  GHZ  state  composed 
of  a  large  number  of  spins. 

We  can  also  directly  calculate  the  metrological  usefulness  of  a  GHZ  state  prepared  in 
the  presence  of  spontaneous  relaxation.  The  favorable  sensitivity  of  the  state  p*  =  p(t*) 
to  rotations  by  angle  fl  around  the  x-axis  can  be  understood  as  the  strong  dependence 
of  the  expectation  value  of  the  parity  operator  fr  = 

P(ft)  =  tr[p*(Q)7r]  (73) 


=  tr  p*  (ft/  cos  Q  —  crj  sin  fl) 


on  the  angle  fl  (where  p*(f2)  results  from  rotating  p*  about  the  x-axis  by  angle  11)  [60]. 
The  phase  sensitivity  of  a  GHZ  state,  denoted  ./#,  is  given  (see  Ref.  [62])  by 

|SP(fi)/an|  I  ..  V'  I  .r.^»TT 


*  = 1  (75) 

V  )  0  j 

where  AP(ll)  =  1  —  P(f2)2  (taking  into  consideration  that  tt'2  =  1)  is  the  uncertainty 
of  the  operator  ft  calculated  in  the  state  p*.  The  approximation  in  Eq.  (75)  is  simply 
that  P(0)  ~  0  and  therefore  AP(0)  ~  1.  This  can  be  checked  explicitly  by  looking  at 
the  large  J\f  limit  of 

P(0)=tr  [?•»]=  (AFT— ■  (76) 


In  the  absence  of  decoherence,  ./#  =  J\f  and  the  Heisenberg  limit  of  phase  sensitivity 
is  obtained.  By  generalizing  Eq.  (58)  to  the  case  where  <r/  is  inserted  on  Af  —  1  of  the 
sites,  we  find  that  in  the  presence  of  decoherence 

=  A fe~Tt/2  Im  [T(2  J  -  2 i7,  t *)Ar"1j  .  (77) 

This  result  is  plotted  in  Fig.  11  for  different  types  of  decoherence.  The  enhancement 
in  ./#  survives  Tx  decoherence  only  if  AfTrt*  <  1  is  satisfied  (the  scaling  by  Af  is 
shown  in  the  inset).  This  result  should  be  contrasted  with  the  effect  of  X2  decoherence 
({Tel,  rud,  rdu}  =  {r,0,0}).  In  this  case  T(2  J,t*)  =  1,  yielding 

=  Afe~rt/2, 


(78) 


CONTENTS 


24 


MTt*  r t* 


Figure  11.  Metrological  gain  over  the  standard  quantum  limit  (..<# /CN)  of  N  =  100 
spins  evolved  under  the  OAT  Hamiltonian  to  a  time  t*  (where  a  GHZ  state  exists  in 
the  absence  of  decoherence).  We  plot  as  a  function  of  decoherence  rates  for  pure 
dephasing  ({Tei,  Tud,  Tdu}  =  {T,  0,0},  red  dotted  line)  and  for  spontaneous  relaxation 
({rel,rud,rdu}  =  {0,r,  0},  blue  dashed  line).  Note  the  two  curves  are  produced  by 
rescaling  the  decoherence  rates  in  different  ways  (in  order  to  show  them  in  the  same 
plot);  if  the  scaling  were  the  same  for  both  plots,  .Of.  would  decay  much  more  quickly 
for  spontaneous  relaxation  than  for  dephasing.  Inset:  Log- log  plot  of  T*,  defined  to 
be  the  decoherence  rate  for  which  Nl  has  decreased  to  N/e,  as  a  function  of  N  (blue 
circles).  The  green  dot-dashed  line  represents  (up  to  a  multiplicative  constant)  1/N 
scaling. 


and  hence  the  precision  enhancement  decays  on  a  timescale  that  is  independent  of  A f 
(consistent  with  results  in  Ref.  [62]).  In  contrast,  the  entanglement  witness  C  decays  at 
an  TV-enhanced  rate  for  either  type  of  decoherence. 

4-3.  Removing  the  effects  of  decoherence  via  measurement-base  feedback 

In  Sec.  4.2,  we  showed  that  spontaneous  relaxation  significantly  degrades  the  precision 
enhancement  of  a  GHZ  state  unless  T  <  J/Af  is  satisfied.  In  this  section,  we  will 
show  that  a  time-resolved  record  of  spontaneous  relaxation  events  provides  sufficient 
information  to  restore  the  phase  enhancement  under  the  much  less  stringent  constraint 
T  <  J.  In  particular,  we  are  imagining  a  situation  where  spontaneous  spin  flips  are 
accompanied  by  the  real  spontaneous  emission  of  a  photon,  such  that  they  can  be 
measured  by  photodetection. 

For  reasons  that  will  become  clear  in  what  follows,  we  take  our  initial  state  to  have 
all  spins  pointing  at  an  arbitrary  angle  9  (rather  than  9  =  n/2,  as  assumed  in  Sec.  4.2). 
Our  goal  is  to  evaluate  the  expectation  value  of  the  operator  <rj  \ of  at  time  t*, 
which  was  accomplished  above  by  appealing  directly  to  Eqs.  (59)  and  (56).  Pursuing 
a  strategy  similar  to  that  employed  in  Sec.  3.5,  we  first  observe  that  n  spins  initialized 
at  9  =  7t/2,  evolving  in  the  absence  of  decoherence  but  in  the  presence  of  a  longitudinal 
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field  of  strength  h,  would  yield  the  phase-sensitivity  enhancement  ft 

h )  =  n  |lm  [e-2iW*(isin2.7t*)n“1]  | 

=  |sin(2/rt*  +  ir (n  —  1)/2)|  . 


(79) 

(80) 


Now  we  calculate  .M  in  the  presence  of  decoherence,  but  we  hold  off  evaluating  the 
multiple  sums  and  integrals  that  yield  the  functions  \fb(s/,t)  in  closed  form,  and  instead 
obtain  .M  by  working  directly  with  Eq.  (48) 

■Jt  =  ^  j  V("R,,  k,  r)Q(s,  7 Z.  r).  (81) 

First,  we  evaluate  Q(s,  R,t),  obtaining 


g{s,  n,  r)  =  (J\f-K)  cos27?  ^0  Im  [e~2ijTg-(s  =  2  Jt  -  2 iyt)^”7*”1] ,  (82) 

which  only  depends  on  the  TZj  and  r3  only  through  their  sums  7 Z  =  J2jTZj  and 
h  =  (J/t)  Tj.  With  the  judicious  choice  6  =  tt  —  2  tan”1  (e27**),  we  can  rewrite 


g  (2 Jt  —  2 iyt) 


i  sin(2  Jt) 
cosh  2yt 


(83) 


and  thus  we  obtain 

*  =  (84) 

The  initial  value  of  9 ,  which  places  the  initial  spins  slightly  above  the  xy  plane,  was 
carefully  chosen  so  that  g~(2Jt  —  2 iy)  oc  sin(2Jt).  This  finite  tipping  angle  is  required 
to  precisely  cancel  the  null  measurement  effect,  which  causes  the  ^-projection  of  each 
spin  to  change  in  time  due  to  the  lack  of  emission  of  a  photon,  and  thus  ensures  that  at 
time  t*  the  unflipped  spins  are  brought  down  into  the  xy  plane.  Finally,  we  can  write 


^  =  ^2  dh  P{K,  h)Jt(N  -  n,h  =  Jr/t),  (85) 

■R  d 

where  P{JZ,  h )  is  obtained  by  carrying  out  in  Eq.  (84)  with  7Z  and  r  =  ht/  J  held 
fixed.  As  in  section  3.5,  we  can  interpret  P( 77,  h )  as  the  probability  to  have  1Z  flipped 
spins  contributing  an  effective  magnetic  field  h  to  the  dynamics  of  the  remaining  spins. 
For  each  particular  value  of  7 Z  and  h,  the  function  (■ n ,  h )  (where  n  —  J\f  —  TZ)  yields 
the  precision  enhancement  of  a  GHZ  state  produced  from  n  spins  in  a  longitudinal  field 
of  strength  h. 

If  an  experiment  can  record  the  times  (t\, . . .  ,tn )  at  which  photons  are  emitted, 
thus  gaining  access  to  both  TZ  and 

j  ^ 

h  =  -  Jj(2 tj  -  t),  (86) 

3= 1 


then  that  experiment  produces  the  conditional  density  matrix  p(n,  h ),  corresponding  to 
a  GHZ  state  of  n  spins  produced  by  one-axis  twisting  in  the  presence  of  the  longitudinal 


ff  Note  that,  for  ft  =  0  and  n  odd,  (n,  ft)  =  0.  This  is  because  the  GHZ  state  created  when  n  is  odd 
is  rotated  from  the  GHZ  states  we  consider  by  an  angle  of  7t/2  in  the  xy  plane. 
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field  h.  The  effect  of  this  field  is  simply  to  rotate  the  system  around  the  2- axis  by  a 
(shot-to-shot  random)  angle  5  =  (ht*+n(n— 1)/2).  However,  rotation  by  a  random  angle 
only  causes  decoherence  if  the  value  of  that  angle  is  not  known.  Because  the  experiment 
measures  8  indirectly  via  the  photon  emission  record,  it  is  possible  to  remove  the  effect 
of  h  in  any  experimental  shot  by  applying  the  rotation  operator  R(8)  =  exp  (— iSz8 )  to 
the  conditional  density  matrix  p(n,h).  In  this  way  we  create 

p(n,  0)  =  R{8)p{n,  h)T$(8),  (87) 

which  is  GHZ  state  of  the  form  in  Eq.  (70)  containing  n  spins,  and  thus  obtain  a  precision 
enhancement  of  n.  Because  the  expected  value  of  n  will  decay  only  at  the  bare  rate  T, 
there  is  no  longer  an  enhancement  by  J\f  in  the  decay  of  precision.  This  measurement- 
based  coherent  feedback  can  also  be  applied  in  the  context  of  spin-squeezing,  where  once 
again  it  can  vastly  improve  the  metrological  usefulness  of  a  state  generated  by  one-axis 
twisting  in  the  presence  of  spontaneous  relaxation. 

Before  concluding,  we  note  that  this  feedback  strategy  could,  in  principle,  be  applied 
to  situations  where  both  spontaneous  relaxation  and  spontaneous  excitation  are  present, 
and  even  when  the  coupling  constants  Jij  are  not  uniform.  However,  in  the  former  case 
it  is  necessary  to  independently  record  the  photon  emission  record  corresponding  to 
excitation  and  relaxation  processes,  and  in  both  cases  it  is  necessary  to  obtain  site- 
resolved  (in  addition  to  time-resolved)  information  about  the  photon  emissions. 

5.  Conclusions  and  future  directions 

In  this  paper  we  have  presented  a  comprehensive  theoretical  toolbox  for  understanding 
far-from  equilibrium  dynamics  in  Ising  models  both  with  and  without  decoherence.  The 
underlying  objects  of  interest  are  unequal-time  correlation  functions,  which  are  then 
used  to  compute  spin  squeezing,  dynamical  response  functions,  entanglement  witnesses, 
and  the  effects  of  dephasing,  spontaneous  excitation,  and  spontaneous  relaxation  on 
the  system  dynamics.  We  believe  these  tools  will  be  of  fundamental  importance  in 
understanding  and  optimizing  a  diverse  array  of  systems  in  which  entanglement  is 
engineered  by  Ising  interactions.  In  particular,  these  tools  enable  the  quantification 
of  detrimental  effects  due  to  system-environment  coupling,  even  when  the  coupling  does 
not  commute  with  the  Ising  interactions. 

The  ability  to  compute  dynamics  in  any  dimension  and  in  the  presence  of  non¬ 
commuting  noise  is  a  particularly  surprising  result;  it  is  well  known  that  the  inclusion  of 
non-commuting  but  coherent  linear  couplings  admits  solutions  only  in  highly  specialized 
geometries,  such  as  ID  nearest  neighbor  chains.  The  key  structures  that  allow 
our  solution  for  the  open  system  to  proceed  are  (1)  the  statistical  independence  of 
decoherence  processes  on  different  sites  and  (2)  a  symmetry  between  the  forward  and 
backward  time  evolution  along  a  closed-time  path,  i.e.  the  Markov  approximation.  It 
would  be  interesting  to  understand  to  what  extent  this  simplification  generalizes  to 
other  models  where  the  incorporation  of  decoherence  would — at  first  sight — appear  to 
be  intractable. 
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Appendix  A. 


The  main  goal  in  this  Appendix  is  to  explicitly  cast  the  series  expansion  for  arbitrary 
observables  in  terms  of  the  time-ordered  correlation  functions  encountered  in  Sec.  2  of 
the  text,  thus  bridging  the  gap  between  Eqs.  (45)  and  (48)  of  the  text.  The  summation 
of  the  series  is  carried  out  later  in  Appendix  B. 

Our  starting  point  is  the  series  expansion  for  the  time-evolution  superoperator 


W(t)  = 


X 


(Mr 


dh%(t  -  tn)M%(tn  -  fn_0  . . .  %(t2  -  b(*i 


(A.l) 


This  leads  immediately  to  the  expression  for  A(t)  given  in  the  manuscript  [Eq.  (45), 
reproduced  here  for  convenience]: 


A{t)  = 
Y,  [  dtr 


r»t2 


dti  tr 


A%(t  -  tn)M%(tn  -  tn- 1) . . .  %{t2  -  H)^0(H)p(0)  .  (A. 2) 


In  order  to  simplify  notation  in  the  following  equations,  we  define  time  dependent  jump 
operators  in  the  Heisenberg  picture  of  the  effective  Hamiltonian 

J(j,a,t)  =  (A.3) 

where  7+(_)  =  rdu(ud).  Defining 

/  °°  _  ft  ft.2 

= y,  y,  y  /  dtn...  (a.4) 

n=0  0  0 

we  can  now  express  A{t)  as 

At)  = 

yjtr  Ae~ltn^J (jn,  an,  tn) . . .  J (A,  «i,  ti)p(0)j7t(ji,  au  h) . . .  j\jn ,  an,  tn)eim^« 

=  S/(  57^ (ji)  t\  )  .  .  .  j\jni  Q"ni  tn)A(t)J' ( jni  ®"ni  tn)  ■  •  •  N (jl;  H)^  •  (A. 5) 


In  the  above  expression,  the  time  dependence  of  the  operator  A  is  defined  as 


A{t)  =  eltn^Ae 


(A.6) 


which  is  distinct  from  the  time  dependence  assigned  to  the  operators  <r“  in  defining 
the  J  (j,  a,  t)  (of  course  this  distinction  vanishes  when  considering  time  evolution 
under  a  Hcrmitian  Hamiltonian).  In  the  final  line  the  trace  has  been  removed  (upon 
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rearrangement,  it  becomes  a  completeness  identity),  and  the  expectation  value  is  in  the 
initial  pure  state  |^(0)). 

The  correlation  functions  in  Eq.  (A. 5)  are  explicitly  closed-time  path  ordered,  with 
the  time  ordering  enforced  in  the  limits  of  integration.  The  remaining  task  is  to  explicitly 
separate  all  time-dependence  of  the  correlation  functions  in  Eq.  (A. 5)  due  to  the  anti- 
Hermitian  part  of  TLcfi,  so  that  we  can  directly  employ  the  results  from  Sec.  2.  Because 
the  anti-Hermitian  part  of  the  effective  Hamiltonian  commutes  with  'H,  this  is  relatively 
straightforward.  As  in  the  manuscript,  we  restrict  ourselves  at  this  point  to  considering 
operators  A  that  can  be  written  as  products  of  spin-lowering  and  spin-raising  operators 
abjJ  ( bj  =  ±1)  on  sites  j  e  rj.  We  can  then  write 


A(t)  =  exp 


JSjTrt 


exp 


-2 


Ait), 


with  A(t)  evolving  in  the  Heisenberg  picture  of  Ti  alone 
A(t)  =  eimAe~im. 


(A.7) 


(A.8) 


Similarly,  we  can  rewrite 

=  (A.9) 

with 


a](t)  =  emta?e-im  (A.  10) 

evolving  in  the  Heisenberg  picture  of  TL. 

Replacing  the  operators  in  Eqs.  (A.7)  and  (A.9)  back  into  Eq.  (A. 5),  we  are 
essentially  ready  to  read  off  the  results  of  a  given  term  from  the  expressions  for 
correlation  functions  in  Sec.  2.  However,  anticipating  that  the  terms  in  the  series  will 
be  site-factorizable,  we  pause  here  to  introduce  some  notation  describing  the  occurrence 
of  jump  operators  belonging  to  a  particular  site.  Because  jump  operators  always  occur 
at  the  same  time  along  the  forward  and  backward  evolution  (and  each  operator  on  the 
backward  path  is  the  complex  conjugate  of  a  corresponding  operator  on  the  forward 
path),  we  only  need  to  describe  the  jump  operators  applied  during  the  forward  time 
evolution.  First,  we  define  TTf  to  be  the  total  number  of  jump  operators  of  type  af 
applied  to  site  j  along  the  forward  time  evolution,  and  IZj  =  TZl-  +  TZj.  We  also  define 
{t{, . . .  ,tJiz  }  to  be  the  set  of  times  at  which  those  jump  operators  are  applied  to  the  site 
j,  and 

[  crj(t)  (A.ll) 

Jo 

4(-i)Rj-”j  (A. 

to  be  the  total  amount  of  time  the  j’th  spin  has  spent  pointing  up  minus  the  time 
it  has  spend  pointing  down,  again  during  the  forward  evolution  only  (note  that 


—  (1  oij)nj  t  2 


IZj 

E 

n—  1 
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fcaj(t )  =  0  Vj  7^  77).  Finally  we  use  the  bold  symbols  7 to  represent  vectors 
of  the  quantities  IZj,  Kj ,  and  Tj,  respectively. 

Now  we  can  write  A(t)  as 


A(t)  = 


u  ■ 


l(tn)A(t)a%(tn) . . .  (ti)  exp 


2  yt 


2 L  aia3 


j£v 


b  (A. 13) 


where 


V(H,k,t)=  \\V3(Nv  kvt2  (A.  14) 

3 

=  J]  (e~v^2  (rud)^“  (Tdu)^+  e-2^)  .  (A.  15) 

3 

The  exponential  of  operators  <jJ  Eq.  (A.  13)  leads  to  the  so-called  null  measurement 
state  reduction:  it  causes  the  jth  spin  to  drift  out  of  the  xy  plane  in  the  event  that  it 
has  not  spontaneously  relaxed  (ay  =  1). 

The  expectation  value  in  Eq.  (A.  13)  is  now  precisely  of  the  form  given  in  Eq.  (14), 
with  the  exception  that  we  must  map  yy  — >  Sjt  =  yy  —  2 i'yt  in  order  to  account  for  the 
term  in  square  brackets.  Thus  we  obtain 


G(s,  n ,  t)  =  (ahai  (h)...  ajnan  (tn)A(t)a%  ( tn ) . . .  d"1  (ti)  exp 


2T«E 

j$V 


QtjCTj 


n  [Pi  +  9j{sjt)]  , 


(A. 16) 
(A.  17) 


3 

where  s  is  a  vector  of  quantities  Sj  =  ipj/t  —  2 iy.  Careful  bookkeeping  reveals  that 


Tj 

2 1  ^  ^  bkJjki 

(A.  18) 

kerj 

■d(r) 

= 

j 

(A.  19) 

0j(Tj) 

=  1 0 

\  ^Tj  Thket)  bkJjk 

if  j  e  rj] 
if  j  £  V, 

(A. 20) 

which  allows  us  to  factorize 

g{8,'R,,T)=  YlGjisj^jNj)  (A.21) 

3 

=  (A. 22) 


3 

thus  completing  the  derivation  of  Eqs.  (48-50)  in  the  manuscript. 
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Appendix  B. 


In  Sec.  3  we  encountered  the  functions 

<I>/ =  J2j.  VJ ’  Ki >  tj ) Gj (sJ’Kj’Tj)  (B • 1 ) 

^j(sjit)  =  "y  j  ^3  ('R'J  >  K‘j  >  Tj  )  ^(1  —  aj)Kj  A  Oij  -  ^  Gj(Sj,lZj,  Tj),  (B.2) 

which  we  now  show  how  to  evaluate.  When  j  G  r],  only  the  7 Zj  =  0  term  in 
survives,  and  we  obtain 

=  e~Vlt/2Pj-  (B.3) 

The  function  is  only  defined  for  j  ^  rj,  so  this  case  does  not  apply  to  it.  We 

next  consider  the  case  j  ^  rj,  and  drop  the  site  index  j  (the  derivation  does  not  depend 
on  the  specific  site  j,  though  the  answer  will  depend  on  s,  which  can  be  reindexed  at 
the  end).  We  begin  with  <h(s,f),  which  can  be  simplified  as 


<h(s,  t)  —  5k,  o  + 


=  e~Frt/2g+(s 


00  rt  rt‘2  \ 

EE  dt-n  ■  ■  ■  /  dt i  V(1 Z,  k,  t)G(s ,  7 Z,  r) 

?=±17?=lV°  / 

[st)  +  ^2^2  [  dtK  ■■■  [  dt{P{lZ,  n,  r)pe 
^  -t~>  J  0  o  0 


(B-4) 


(B.5) 


(rud)^“ (rdu)^+  p  f  dtn...  r dt  1  eisr  (B.6) 

'TJ  J  0  io 


=  e~r’,/2g+{st)  +  e_r'"'2  ^  X(i 


(B.7) 


We  will  carry  out  the  integrals  first,  for  which  it  is  convenient  to  treat  the  cases  1Z  even 
and  7 Z  odd  separately.  Remembering  that 


r  =  (1  —  a)  /  az(t ) 

Jo 

—  n  [t  —  2 tK  +  2t7^_i . . .  +  (— l)^2ii]  , 
and  defining  a  new  index  /i  satisfying 

f  S=1  if  -R  is  odd 


(B-8) 


(B.9) 


jf  7^  js  even, 


we  obtain 


*(k)  = 


t  (1  —  k  cos  6) 


(rr  -4^7  + 


t-Kidz\  ( rtV1  j^zt) 


/i  +  l  J  \2^  /  /i! 


(B.10) 


Here  the  j/t  are  spherical  Bessel  functions,  and  the  integral  has  been  carried  out  by 
changing  variables  to  allow  evaluation  of  all  but  one  integral  while  r  is  held  fixed, 
and  then  evaluating  the  remaining  integral  over  r  (this  is  where  the  Bessel  functions 
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arise).  The  remaining  sum  can  be  evaluated  using  relations  obtained  from  the  generating 
function  for  spherical  Bessel  functions 


Si  (x,y)  =  Yx' 


My) 

li=0  ^  +  1)] 

=  sine  a Jy1  —  2 xy 


and 


^2  (x,y)  =  Y 


=  2^x 

fi= 0 


ojM 

(TV- 


cosy 


=  —  cos  \/y2  —  2xy 

xy  xy 

In  terms  of  these,  we  finally  obtain 


X(k)  = 


t  (1  —  K  COS  6) 


(  rt 


V  2s 


(Tr  -  4^7)  Sx  — ,  st  +r(t-  indz)  S2(  —,zt 


rt 


2z ' 


(B.  11) 
(B- 12) 

(B.13) 
(B- 14) 

(B.  15) 


Tedious  but  straightforward  algebra  leads  to  Eq.  (56)  of  Sec.  3.  Notice  that  we  can 
just  as  easily  evaluate  the  similar  expression 


T(s,f)  = 


/  _ ^  °°'  ft  ft2  \ 

dn,o  +  IE  ( tt-ji  ■  ■  ■  /  dti  )  V(1Z,  k, 

\  k=±1TZ=1^°  / 


r)  (  (1  -  a)K  +  a1-—  )  Q(s,  TZ,  r) 


—  e  ritOn 


9  (s 


rt  ft2 

'i  +  EE  /  dt-Ti . . .  /  dt{P(Tl,  k,  T)npe 
Ti  Jo  J  0 


K  1Z 


—  e  —'-g  [St)  +  e  Frt/2^T( 


=  e  Tlt/2g 

giving  rise  to  Eq.  (59)  of  Sec.  3. 


K)K 


(B.16) 

(B.17) 
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